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We use Genetic Algorithms to extract information from several cosmological probes, such as the 
type la supernovae (Snla), the Baryon Acoustic Oscillations (BAO) and the growth rate of matter 
perturbations. This is done by implementing a model independent and bias-free reconstruction of 
the various scales and distances that characterize the data, like the luminosity cLl{z) and the angular 
diameter distance c?a(z) in the Snla and BAO data, respectively, or the dependence with redshift 
of the matter density Q m (a) in the growth rate data, fa$(z). These quantities can then be used to 
reconstruct the expansion history of the Universe, and the resulting Dark Energy (DE) equation of 
state w(z) in the context of FRW models, or the mass radial function f?Af(r) in LTB models. In 
this way, the reconstruction is completely independent of our prior bias. Furthermore, we use this 
method to test the Etherington relation, ie the well-known relation between the luminosity and the 

£N| angular diameter distance, r\ = jjj^z^jyj i which is equal to 1 in metric theories of gravity. We 

find that the present data seem to suggest a 3-cr deviation from one at redshifts z ~ 0.5. Finally, 
we present a novel way, within the Genetic Algorithm paradigm, to analytically estimate the errors 

h^- on the reconstructed quantities by calculating a Path Integral over all possible functions that may 

contribute to the likelihood. We show that this can be done regardless of the data being correlated 
or uncorrelated with each other and we also explicitly demonstrate that our approach is in good 

£SJ agreement with other error estimation techniques like the Fisher Matrix approach and the Bootstrap 

Monte Carlo. 
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A decade after the fundamental discovery of the acceleration of the Universe [T] , cosmologists are still puzzled as to 
the nature of the dark energy supposed to be responsible for such acceleration. Many alternatives have been proposed: 
a cosmological constant A (or more generally an approximately constant vacuum energy); an effective scalar field (like 
in quintessence [2], [3] etc.); various relativistic ether fluids (like e.g. Chaplygin gas[4]); several ad hoc modifications 
of general relativity on UV scales (for instance scalar-tensor theories [5], f{R) gravity[6], Gauss-Bonnet terms[7], 
^NJ Torsion cosmology[H], etc.); extra symmetries like conformal Weyl gravitv[5].[TU] or massive gravitons[dl ; introduction 
of extra dimensions (like in DGP models 12J, or Kaluza-Klein compactifications) ; new effective interactions (like 
the Galileon[T3] field, etc.); finally there is the possibility that the observed dimming of supernovae is not due to 
acceleration in a FRW universe but to large inhomogeneities from gravitational collapse, or from large LTB Voids [TH- 
IT7] . For reviews on many of these interesting alternatives we refer the interested reader to Ref. [IB] , 

All these alternatives could in principle be distinguishable if we had detailed knowledge about the redshift evolution 
of a handful of parameters associated with concrete physical quantities, like the matter content fl m (a), the rate of 
expansion H(a), the luminosity distance d,L(a), the angular diameter distance d,A(a), the galaxy and cluster number 
counts dN/dQ(a), the deceleration parameter q(a), the cosmic shear S(a), the density contrast growth function /(a), 
• • and 7(a), the Jeans length of perturbations c 2 s {a), the anisotropic stresses of matter 77(a), the scalar torsion parameter 
^ s(a), the bulk viscosity ((a), e te- 

Fortunately, we now live in an era where all these parameters could in principle be measurable in the not too 
far future thanks to the proposed probes of LSS and the CMB, that will determine the following observables with 
increasing accuracy: the matter power spectrum P(k,a), the Snla distance modulus ^t(a), the Baryon Acoustic 
Oscillation scale, #BAo( a )> the cluster number counts dN c /dVl(a), the galaxy and halo mass functions n(< M), the 
lensing magnification and convergence /i, k, the Redshift Space Distortions and bias /3(a), b(z) of LSS, the local rate of 
expansion Ho, the age of the Universe to, the CMB temperature and polarization anisotropics Ci(TT,TE, EE, BB), 
the Integrated Sachs- Wolfe and the Sunyaev-Zeldovich effects, and possibly also the fractal dimension of space time. 

There is, however, a limit to what we can say about the physics responsible for acceleration from observations. How 
many parameters can we constrain? What is the optimal parametrization of the linear perturbation equations? How 
much can we extract from nonlinear regime? Can we interpolate between super-horizon scales sub-horizon mildly- 
nonlinear and full nonlinear scales? Can we parameterize wide classes of models? What is the role of systematics on 
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uncertainties? All these are crucial questions that will have to addressed in the near future when we have a significant 
increase in data from LSS and CMB. For the time being, we should be as open as possible to any of the above 
alternatives. In that respect we would like to pose the appropriate questions to the right number of cosmological 
observables without having to assume a model, i.e. without any prejudices as to the underlying space-time structure 
and dynamics. 

In particular, we don't want to determine the likelihood contours for the N-dim parameter space of a fully specified 
ACDM model, but rather approach the observational data sets without prejudices, in a formal way, and then use the 
outcome to constrain whatever model or realization we may envision, be it mCDM, or DGP, or f(R), or LTB, etc. 
This is the reason we develop the Genetic Algorithm (GA) approach to DE modeling. Genetic algorithms have been 
used in the past in Cosmology, see [HI [20] , to address supernovae distance modulus fitting in a model independent 
way, but this is the first time it is done in full scale to the multiprobe space of Dark Energy modeling. 

In Section [TT] we describe the Genetic Algorithms as applied to the cosmological problem at hand and we also 
describe a novel prescription for the computation of errors in the context of GA by doing a path integral over all 
possible functions that may contribute to the likelihood. In Section [Til] we give an account of the data sets used in this 
paper, the type la supernovae (Snla) data, the Baryon Acoustic Oscillations (BAO) data and the growth rate data, 
used for the computation of the redshift dependence of a handful of important parameters like the luminosity and 
angular diameter distances <Il{z) and cIa(z) respectively and the matter density f2 m (z). In Section IV we describe 
the results of the GA analysis in the context of various models: FRW/DE, LTB, and modified gravity theories. In the 
context of FRW/DE models, we compare many different diagnostics like the equation of state w(z), the deceleration 
parameter q(z) and the Om statistic and we find that the best in constraining the evolution and properties of DE 
is the deceleration parameter q(z) as it requires no prior assumptions or a priori set parameters like w{z) does and 
has the smallest errors especially at small redshifts. In the case of the LTB models, we reconstruct the mass radial 
function Ojf(r) and we find that it is in some tension with the commonly used profiles found in the literature. For the 
case of the modified gravity theories we reconstruct the important parameter G c ff , which is a smoking gun signature 
for deviations from GR if it is different from unity at any redshift, and we find a 3er bump at intermediate redshifts. 

V and find a similar deviation in 



We also include a test of the Etherington relation r\ = = 1 in Section 

the same redshift range. However, since the statistical significance is not high enough we cannot draw any strong 
conclusions. Finally, in Section [VI| we give our conclusions. 



II. GENETIC ALGORITHMS 



A. Theory 

In this section we will briefly introduce the Genetic Algorithms (GA). For a more detailed description and the 
application of G As to cosmology we refer the interested reader to Refs. [T5] and [50] ■ The GAs are algorithms loosely 
modeled on the principles of the evolution via natural selection, where a population of individuals evolves over time 
under the influence of operators such as mutation (random change in an individual) and crossover (combination of 
two or more different individuals), see Fig. [I] for a schematic description of the operations. The probability or in 
other terms the "reproductive success" that an individual will produce offspring (the next generation) is proportional 
to the fitness of this individual, where the fitness function measures how accurately it describes the data and in our 
case it is taken to be a x 2 - 

The algorithm begins with an initial population of candidate individuals (in our case functions), which is randomly 
generated based on a predefined grammar of allowed basis functions, eg sin, exp, log etc and operations +, — , x, etc. 
In each successive generation, the fitness functions for the individuals of the population are evaluated and the genetic 
operators of mutation and crossover are applied in order to produce the next generation. This process iterated until 
some termination criteria is reached, e.g. a maximum number of generations. To make the whole process more clear 
we can also summarize the various steps of the algorithm as follows: 

1. Generate an initial random population of functions M(0) based on a predefined grammar. 

2. Calculate the fitness for each individual in the current population M(t). 

3. Create the next generation M(t + 1) by probabilistically choosing individuals from M(t) to produce offspring 
via crossover and mutation, and possibly also keeping a proportion of the previous generation M(t). 

4. Repeat step 2 until a termination goal has been achieved, eg a maximum number of generations. 

We should point out that the initial population M (0) depends solely on the choice of the grammar and the available 
operations and therefore it only affects how fast the GA converges to the best-fit. Using a not optimal grammar may 
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FIG. 1: Left: The crossover operation: the "chromosomes" of the parents are shown by yellow and red circles, where each 
individual circle indicates a gene. Each parent's chromosome is split at the same point and the remaining parts are exchanged. 
Right: The mutation operation: the "chromosome" of a parent is shown by the yellow circles, where each individual circle 



indicates a gene. One specific gene is chosen at random and its value changes (from yellow to pink). See Section II B for more 
details and a hands-on example. Both images from Ref. [21] , 



result in the algorithm not converging fast enough or being trapped in a local minimum. Also, two important 
parameters that affect the convergence speed of the GA are the mutation rate and the selection rate. The selection 
rate is typically of the order of 10% ~ 20% and it affects the number of individuals that will be allowed to produce 
offspring. The mutation rate is usually of the order of 5% ~ 10% and it expresses the probability that an arbitrary 
part of individual will be changed. If either of the two rates is much larger than these values then the GA may not 
converge at all, while if the two rates are much smaller then the GA will converge very slowly at some local minimum. 

The difference between the GAs and the standard analysis of observational data, ie having an a priori defined model 
with a number of free parameters, is that the later method introduces model-choice bias and in general models with 
more parameters tend to give better fits to the data. Also, GAs have a clear advantage to the usual techniques when 
the parameter space is large, far too complex or not well understood, as is the case with dark energy. Finally, our goal 
is to minimize a function, in our case the x 2 > not using an a priori defined model, but through a stochastic process 
based on a GA evolution. In this way, no prior knowledge of a dark energy or modified gravity model is needed and 
our result will be completely parameter-free. 

In other words, the GA does not require us to arbitrarily choose a DE model, but uses the data themselves to find 
this model. Also, it is parameter free as the end result does not have any free parameters that can be changed in order 
to fit the data. So, in this sense this method has far less bias than any of the standard methods for the reconstruction 
of the expansion history of the Universe. This is the main reason for the use of the GAs in this paper. The Genetic 
Algorithms were first implemented in the analysis of cosmological data in Refs. [T(5] and [217] , 

In the present analysis the actual fitting of the data was done with a Mathematica code written by one of the 
authors and a modified version of the GDF v2.0 C++ program 1 developed by I. Tsoulos et al [22] . 



B. An example 

In this Section we will briefly describe a simple example 2 of how the GA determines the best-fit given a dataset. 
We will not describe any of the technicalities but instead we will concentrate on how the generations evolve over 
time, how many and which individuals are chosen for the next generation etc. Also, in order to keep our example 
as simple as possible we will assume that our grammar includes only basic functions like polynomials x, x 2 etc, the 
trigonometric functions sin(x), cos(x), the exponential e x and the logarithm ln(x). As mentioned in the previous 
section, the first step in the GA is setting up a random initial population M(0) which can be any simple combination 
of the grammar and the allowed operations, eg Jga,i(x) — ln(x), fcA.2(x) = —l + x + x 2 and fGA,z{x) = sin(x). The 
number of functions in the genetic population is usually on the order of a few hundreds. 

Next, given some data (xi,yi, a{) the algorithm measures the fitness of each solution by calculating their x 2 defined 



1 Freely available at http://cpc.cs.qub.ac.uk/summaries/ADXC 

2 This is only a schematic description of how the GA works and this is done solely for the sake of explaining the basic mechanisms of the 
GA. 
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as 

^^fy^lM) 2 (2 .d 

i=l ^ 1 ' 

where for our hypothetical example Xi = 200, x| = 500 and x§ = 1000. The selection per se is done by implementing 
the "Tournament selection" method which involves running several "tournaments" , sorting the population with respect 
to the fitness each individual and after that a fixed percentage, adjusted by the selection rate (see the previous 
subsection), of the population is chosen. As mentioned before, the selection rate is of the order of 10% ~ 20% of 
the total population, but lets assume that the two out of the three candidate solutions (/ga.iW and /ga,2(x)) are 
chosen for the sake of simplicity. 

The reproduction of these two solutions will be done by the crossover and mutation operations. The crossover will 
randomly combine parts of the "parent" solutions, for example this may be schematically shown as 

faA,i(x) ® foA,2( x ) -> Uga.i{x)J GAi2 {x)J ga ^{x)) 

= (ln(a; 2 ),-l + ln(a; 2 ),-l + ln(x)) 

In the next step the GA will implement the mutation operation. The probability for mutation is adjusted by the 
mutation rate, which as we mentioned is typically of the order of 5% ~ 10%. In this example this may be a change 
in the exponent of some term or the change in the value of a coefficient. For example, for the candidate solution 
/ga,3{x) this can be schematically shown as Jga,3(x) = — 1 + ln(ir) — > — 1 + ln(a; 3 ), where the power of the x term 
was changed from 1 to 3. 

Finally, at the end of the first round we have the three candidate solutions 

M(l) - (f G A,i(x), Jga^x), f GA , 3 (x)) = 
(\n{x 2 ),-l + \n(x 2 ),-l + \n(x 3 )) 

In the beginning of the next round the fitness of each individual will again be determined, which for our population 
could be (Xi,X2'Xi) = (150,300,400), and the selection and the other operations will proceed as before. It is clear 
that even after just one generation the combined x 2 01 the candidate solutions has decreased dramatically and as we 
will see in the next section, it usually takes a bit more than 100 generations to reach an acceptable x 2 and around 
500 generations to converge on the minimum. After a predetermined number of generations, usually on the order of 
1000, has been reached or some other goal has been met, the GA will finish, with the best-fitting function of the last 
generation being considered as the best-fit to the data. 



C. Error analysis 

1. Uncorrelated data 

One possible shortcoming of the Genetic Algorithms is that they do not directly provide a method to estimate the 
error on the derived best-fit function. One way to around this issue was used in Ref. [2U], were a bootstrap Monte 
Carlo was performed in order to estimate the error of the derived best-fit. We refer the reader to Refs. [2D] and [53] 
for more details on the bootstrap approach to error estimation. In this paper we will follow a more direct approach 
by deriving some analytic results. First we will remind the reader about some important facts regarding the normal 
distribution. If the normal distribution with a zero mean and a a 2 variance is given by 

f(x -° ) = vk^(-&)- (2 - 2) 

Then in the case of one variable only the la confidence interval (CI) or equivalently the region that contains approx- 
imately the 68.3% of the values is found by calculating the area within x G [—If, la] or 

CI(la) = J dxf(x,a)=J^ dx-^=-exp(-^j = erf (1/V2) , (2.3) 

where erf(x) is the error function. This can easily be generalized to n as as CI(na) = erf (n/\/2). In our case, our 
likelihood functional is 



£ = AAexp(- X 2 (/)/2) 



(2.4) 
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where f(x) is the function to be determined by the Genetic Algorithm, x 2 (f) 1S the corresponding chi-squared for 
N data points (xi, yt, <7j), defined Eq. (2.1 1. Then, determining the normalization constant Af is more complicated 
than that of a simple normal distribution, as we have to integrate over all possible functions f(x) or in other words 
perform a "path integral" 

"®f £ = J®f Mcxp(- x 2 (f)/2) = 1, (2.5) 

where 23/ indicates integration over all possible values of f(x). The reason for this is that as the GA is running it may 
consider any possible function no matter how bad a fit it represents, due to the mutation and crossover operators. 
Of course, even though these "bad" fits will in the end be discarded, they definitely contribute in the total likelihood 
and have to be included in the calculation of the error estimates of the best-fit. 

The infinitesimal quantity Df can be written as Df = YiiLi 4f<i where dfi and /, are assumed to mean df(xi) and 
f(xi) respectively, and we will for the time being assume that the function / evaluated at a point Xi is uncorrelated 
(independent) from that at a point Xj. Therefore, Eq. (2.5) can be recast as 



fs)fC= f JJd/iATexpf-^ 

J J -°° s=i V i=i 

N r +co / 1 / f \ 2N 

= I]/ dfrMexpl- 1 ' yi ~ fi 

i=l J-oo \ 
N 

= Af-(2ir) N/2 l[a 



a, 



i=l 



which means that N = (^(2ir) N ^ 2 YiiLi <J i S j ■ Therefore, the likelihood becomes 



C = ~„ ^m^N ex P (-X 2 (f)/2) (2.6) 



1 

(2 7 r) /v/F nr= 1 * 

or if we take into account our assumption that the function / evaluated at each point Xi is independent 



i /i / y%- h 



where 



i / i / 'vi - /i x J 



exp 



(27t) 1/2 ( j 4 " V 2 V o-i 
Now, we can calculate the la error 5fi around the best-fit fbf(x) at a point X4 as 

Jhf{*i)-5h (2^) 1/ V i V 2 ^ CTi 



If we demand that the errors Sfi correspond to the la error of a normal distribution, then from Eq. (2.9 1 we can solve 
the following equation for Sfi numerically 



CI{xi,8fi) = evi (l/VS), (2.10) 

and therefore determine the la error 5fi of the best- fit function fbf{x) at each point X4. However, this will lead to 
knowledge of the error in specific points Xi, which is not ideal for our purpose, ie to have a smooth, continuous and 
diffcrentiable function. Therefore, we will create a new chi-square defined as 

N 2 

Xli{5fi) = (Cj(*i,*/i) - erf (l/V2)) (2.11) 



() 



and we will also parameterize Sf with a second order polynomial Sf(x) 
combined chi-squared x 2 {fbf + Sf) + Xci(Sf) f° r the parameters (a, b, < 



is given by Eq. (|2.11|). 



Then, the la region for the best-fit function fbf(x) 
Sf(x)]. 



a + bx + ex 2 . Finally, we min imize the 
where x 2 is given by Eq. (2.1 1 and Xci 



will be contained within the region 



Sf(x),f b f{x) 

The reason why we perform a combined fit of both chi-squares is that we need to ensure that the derived error 
region corresponds to small perturbations around the best-fit and not to the whole infinite dimensional functional 
space usually scanned by the genetic algorithm. This is inline to what is usually done in traditional error estimates 
either with the Fisher matrix approach, where the errors are estimated from the inverse Fisher matrix, ie the covariance 
matrix, evaluated at the minimum, or in Monte Carlo simulations, where the errors are determined by sampling the 
parameter space around the best- fit |23) . 



2. Correlated data 



Our results can also be generalized even in the case when the data are correlated with each other. In this case the 
chi-square can be written as 



JV 



(ift - /(*<)) cr/ fa - f( Xj )) 



(2.12) 



where Cjj is the covariance matrix. The derivation proceeds as usual but now we will have to rotate the data to a basis 
where they are not correlated with each other. To do so we define a new variable Si = (jjj — f(xj)), where can 
be found by decomposing the inverse covariance matrix C" 1 = D T D by using Cholesky decomposition and then the 
chi-squared can be written as x 2 = s 2 . Then, we have that dsi...dsN = \D\ dfi...dfj^ and the functional integration 
can proceed as usual and the normalization can be found. Unsurprisingly, the resulting normalized likelihood can be 
written as 



C = 



1 



(2tt) 7V/2 |C| 1/2 



exp i 



-X 2 (/)/2) 



(2.13) 



with x (/ ) give n by Eq. (2.12). In the limit wh ere the points become uncorrelated, ie is diagonal with Cu = erf, 
then Eq. (2.13) clearly agrees with that of Eq. (2.6) as in this case | C | 1//2 = Y[f <Ji- 

With this in mind we can also calculate the ler confidence interval (CI) for the best-fit function fbf(x). The 
procedure is similar as in the uncorrelated case, but now we must now rotate to the new basis Sj = (?/j — f(xj)). 
However, now instead of integrating over (— oo, oo) we will instead integrate over the region [fof — Sf, /&/ + Sf] 



CI= QfC = 



hf+Sf N 

ht-'f ti (2tt) 



1 



N/2 ,£,1/2 



ex P (-x 2 (/)/2) 



with x given by Eq. (2.12). Rotating to the new basis Si = (yj — f(xj)) we have 



ds±...dsN 

Sbf,i 

Ssi 
LDI 



\D\dh...df N 

DijiVj- fbf(xj)) 
-DijSfj 

|cr 1/2 



and Eq. (2.14) becomes 



(2.14) 



(2.15) 
(2.16) 
(2.17) 
(2.18) 



CI 



n 



ds. 



1 



-6s, ' (2tt) 1/2 
where, after doing a gaussian integral, the Cli is given by 



exp 



CL = 



erf 
erf 



Ssi + s b f,i 



erf 



Ssi - Sbf, 



V2 



'>^v; n...; • sn 

V2 



N 



licit 



-Dyf-yj+ftrfj+^j) 



(2.19) 



(2.20) 



7 




Real: a+(x-b) Exp[-c x A 2], (a,b,c)=(0.25,0.25,0.5) 
Best-fit (a,b,c)=(0.189,0. 177,0.458), ^,=31.80 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 



Fisher Matrix Method 




FIG. 2: Comparison of the different metho ds for the error estimation. Top left: The real (green line) and best-fit to the mock 
data (red line) of the model of Eq. (2.21 1. Top right: The estimate of the error region by a Fisher Matrix approach. Bottom 
left: The estimate of the error region by a Bootstrap Monte Carlo approach approach (for more details see the text). Bottom 
right: The estimate of the error region by the Path Integral approach. 



For uncorrelated N measurements we have that C%j = diag (erf, ...,a%) and therefore Dy = ±diag . Then 



by choosing the minus sign for Dy, Eq. (2.20) reduces to Eq. (2.9 1 as expected. In order to actually calculate the error 



Sf(x) for the best-fit function fbf(x), we can now follow the same pro cedur e as before b y min imizing the combined 
chi-squared X 2 (fbf + $f) + Xcii^f) now with the \ 2 given by Eq. (2.12) and in Eq. (2.11) the Glj, instead given 
by Eq. plo) ). 



3. Comparison with other methods 



In order to assess how rigorous the path integral method described above is, we will now consider an explicit example 
and compare it numerically with a Bootstrap Monte Carlo and the Fisher Matrix approach. 

First, we created a set of 31 mock data points (ij, j/,, a yi ) with noise in the range x £ [0, 1.55] based on the model 



f(x) — a + (x — 6) exp(— ex 2 ), 



(2.21) 



where the parameters (a, b, c) have the values (0.25, 0.25, 0.5) respectively. The choice of the model was completely ad 
hoc, except for the requirement to be well behaved (smooth) and that it exhibits some interesting features, such as a 
maximum at some point x. This was done in order to test if our methodology is capable of detecting such features 
and if it c an estimate the correct values for the errors at the same time. Then we fitted the model f(x;a,b,c) of 

Eq. (2.21) back to the mock data by minimizing the x 2 (a,fr, c) = ^ ^ v '~f(*" a ' b ' c ) ' s J . The best-fit was found to be 

(a,6,c]~„ = (0.189 ± 0.170, 0.177 ± 0.178, 0.458 ±0.103) with a x 2 m in = 31 - 80 - The original or "real" model along the 
best fit for the same model are shown in the top two plots of Fig. [2j The error region on the top right plot of Fig. 
[2] was estimated by following a Fisher Matrix approach, where the error of the best- fit function / can be calculated 



by <Jf{x) 2 — J2i jCijdif(x;a,b,c)djf(x;a,b,c)\ min evaluated at the best-fit, as described in [53], while the dummy 
variables correspond to the three parameters (a,b,c). The covariance matrix can be simply calculated as the 
inverse of the Fisher matrix Cy = F^ 1 where = \dij\ 2 {a, b, c)\ m in evaluated at the best-fit. 

For the Bootstrap Monte Carlo, we created 100 mock data sets from the original or "real" data with replacement 
and then applied the GA on them to determine the best-fit for each one [23]. After following the prescription as 
described in Ref . [TH] , we calculated the error region and we present it in the bottom left plot of Fig. [2j There seems 
to be a sweet spot in the error region at x ~ 1, ie very small errors, but this is due to the stochastic nature of the 
Bootstrap and the behavior (sharp transition) of the data at this region, ie "below" the best-fit for x < 0.9 and above 
the "best-fit" at x > 0.9. 

Finally, we also calculated the error region based on the Path Integral approach described in the previous sections. 
In this case, the best-fit by the GA was found to be 

fGA{x)= X x , (2.22) 
X + 

for Xmin = 30.76, which clearly is much better than the best-fit of the real model itself! This obviously reinforces our 
belief that the GAs can provide really good descriptions of the underlying data. However, it should be noted that 
extrapolating the best-fit GA function outside x S [0, 1.55], ie in the region x > 1.55, is not advised since the lack of 
data in that region will lead to dramatic deviations of the GA prediction and the real model. 

The best-fit f(x) function and its la errors are shown in the bottom right plot of Fig. [2] As it can be seen, 
the agreement between all three methods is remarkably good, something which lends support to our newly proposed 
method for the error estimation in the GA paradigm. Furthermore, we should note that since only the Bootstrap 
Monte Carlo is applicable to the case of the GAs (as the best-fit has no parameters over which one may differentiate, 
like one would do with the Fisher Matrix), our method has the obvious advantage that it is much faster that the 
Bootstrap, less computationally expensive and equally reliable. 



III. THE DATA 



A. Supernovae of type la 

The analysis of SN la standard candles is based on the method described in Ref. [21]. We will mainly use the 
recently released update to the Union 2 set, i.e. the Union v2.1 dataset [25] . 

The SN la observations use light curve fitters to provide the apparent magnitude m(z) of the supernovae at peak 
brightness. This is related with the luminosity distance di,(z) through m(z) = M + 51og 10 (c?i/10pc), where M is the 
absolute magnitude. In flat FRW models, the luminosity distance can be written as di,{z) = (l + z ) f h(x)/h ' wnere 
H(z) = | is the Hubble parameter and a{t) is the scale factor of the FRW metric. In LTB models, the luminosity 
distance can be written as dr,{z) = (1 + z) 2 A(r(z), t(z)), where A(r, t) is the scale factor but which now also depends 
on the radial coordinate. In the FRW limit we have A(r, t) = r ■ a(t). 



TABLE I: The datasets used in this analysis and the information that can extracted by using the Genetic Algorithms. In the 
last two columns we show to what exactly this information corresponds to in generic FRW/Dark Energy and LTB models. For 
more details see also the text. 

Dataset Information FRW/DE models LTB models GA function 

Snla cIl(z) 
BAO H{z),d A {z) 
Growth-rate / = 7,flm(a) 

"where in both cases D v (z) = ((1 + z) 2 D A (z) 2 y^j) ^ and H{z) = H L (z) for the LTB 
'where fl^f (r) is the mass radial function. 



(l + ^)/o Z j|) (l + z) 2 A(r(z),t(z)) 

(l + z)£(z)^giWfi) a G ^ 2 ( 2 ) 
b 



GAi{z) = H d L (z) 

zd z (z) 



'3^0(^0111=1) 

GA 3 (z) =/a8(z) 
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TABLE II: The BAO data used in this analysis. The first six data points are volume averaged and correspond to Table 3 of 
[30| . Their inverse covariance Matrix is given by (3.8 1. 





6dF 


SDSS 


WiggleZ 


z 


0.106 


0.2 0.35 


0.44 0.6 0.73 


d z 


0.336 


0.1905 0.1097 


0.0916 0.0726 0.0592 


Ad z 


0.015 


0.0061 0.0036 


0.0071 0.0034 0.0032 



Defining the dimensionless luminosity distance as Dl(z) = cIl(z)/H , the theoretical value of the apparent 
magnitude is 

rrh h (z) = M(M,H )+5log w (D L (z)), (3.1) 

where M = M - 5 log 10 h + 42.38 [53] . 

The theoretical model parameters are determined by minimizing the quantity 



XsN] 



i=l 



u^^ r^t.^ ) > (^) 



where N is the number of the SN la dataset, pj is the set of parameters to be fitted, and i are the errors due to 
flux uncertainties, intrinsic dispersion of SN la absolute magnitude and peculiar velocity dispersion. The theoretical 
distance modulus is defined as 

lhb{z%) = ™>tb{zi) - M = 51og 10 (D L (z)) + fi Q , (3.3) 



where /irj = 42.38 — 5 log 10 h. More details for the usual minimization of (3.2 1 are described in detail in Refs. 
Therefore, by fitting the Snla data with the Genetic Algorithm we directly reconstruct the luminosity distance D^, 
which from now on we will call GA±(z) = Dl(z). At this point it should be noted that in determining the best-fit we 
don't make any a priori assumption about the curvature of the universe or any dark energy model, but instead we 
allow the GA to determine the luminosity distance Dz, directly from the data and the result is shown in Fig [3] As it 
can be seen, the prediction of the GA is very close to best-fit flat ACDM model with fi m = 0.27. 



B. Baryon Acoustic Oscillations 

The BAO data are usually given, in FRW models, in terms of the parameter d z {z) = 1ba ^^ 9 ^ , where lBAo{zdrag) 
is the BAO scale at the drag redshift, assumed known from CMB measurements, and [35J 

D v {z)=U + zfD A {zf^j' 3 (3.4) 
is the usual volume distance. From this it is easy to see that for z <C 1 d z scales as d z ~ 1bao . 



where 



In LTB models one has to take into account the effect of the inhomogeneous metric, which induces and extra factor 
(1 + z)^(z). Therefore, for these models [H] 

d z (z) = (l + z)Z(z) — — , (3.5) 

Uv{z) 

= ( A'(r(z),t(z)) A'( roo ,t e ) \ 1/3 ( A(r(z),t(z)) A^Q \ 2/3 

?l > \ A'(r(z),t e ) A'( roo ,to)J V A{r{z),t e ) A( roo ,t ) J ' 1 ' ' 

using a suitable early time t e = t(z ~ 100). 

In this analysis we use the 6dF, the SDSS and WiggleZ BAO data shown in Table 2. The WiggleZ collaboration 
[30j has measured the baryon acoustic scale at three different redshifts, complementing previous data at lower redshift 
obtained by SDSS and 6DFGS [33]. 
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FIG. 3: Up: The GA best-fit distance modulus fi compared to the best-fit flat ACDM model (f2 m = 0.27). Down: The evolution 
of the GA best-fit di,(z) with redshift compared to ACDM for different values of fi m . Not surprisingly the prediction by the 
GA is very close to the ACDM best-fit fi m = 0.27. 

The chi-square is given by 

xIao = E( d * - d (*i)) C 7j L Vi - <%))> (3-7) 
where the indices i,j are in growing order in z, as in Table For the first six points, C^ 1 was obtained from the 
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FIG. 4: The GA best-fit d„(z) compared to ACDM (!T2 m = 0.27). 



covariance data in [30] in terms of d z : 



/ 4444 


0. 


0. 


0. 


0. 


0. \ 


0. 


30318 


-17312 


0. 


0. 


0. 


0. 


-17312 


87046 


0. 


0. 


0. 


0. 


0. 


0. 


23857 


-22747 


10586 


0. 


0. 


0. 


-22747 


128729 


-59907 


V o. 


0. 


0. 


10586 


-59907 


125536 / 



(3.8) 



In order to simplify the analysis we will actually use the function 

d z (z) = lBAo{n ° ai = 1) GA 2 (z) (3.9) 

z 

where the function GA 2 (z) is to be determined by the Genetic Algorithm and and the factor 1/z was introduced in 
order to remove the singularity at low z. The reason we use the constant lBAo( z drag) for fiom = 1 is twofold: first, one 
should note that Ibao just shifts d z vertically and the difference with the "real" value can be easily accommodated 
by the GA and second, we also want to use the best-fit to reconstruct the LTB as well, which is actually a matter 
dominated model. In Fig. [4] we show the GA best-fit d z (z) compared to ACDM for fi m = 0.27. 



C. The growth rate of density perturbations 



The growth rate data are given in terms of the combination of parameters fa8(z) = f(z) ■ (7s(z) where f(z) is the 
growth rate of matter perturbations 

f (a) = -rr-f 3.10 
a ma 

and S m = is the linear growth factor which, for models where dark energy has no anisotropic stress, satisfies the 
following ODE 

5 + 2HS = 4nG NPm S Pm (3.11) 
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The exact solution of Eq. (3.11 1 for a flat GR model with a constant dark energy equation of state w is given for the 
growing mode by [31-55]: 



5(a) = o-aFif-i , i - J-; 1 _ A ; a -3» (1 _ 

\ 3w 2 2w; 6u> / 

for H 2 (a)/H 2 = n m a~ 3 + (1 - ft m )a- 3(1+to) , (3.12) 
where 2Fi(a, 6; c; z) is a hypergeometric function defined by the series 

2 F l (a,b;c;z = TV 1 \ I z 3 - 13 

1 (all (o) ^— ' l(c + n n! 

on the disk |z| < 1 and by analytic continuation elsewhere, see Ref. [31] for more details. Also, a$(z) corresponds to 
the rms mass fluctuation and is given by 

a 8 (z) = a 8 (z = 0) 6rn{z) (3.14) 
d m (z = 0) 

From the definitions of f(z) and a$(z) it is easy to prove the following 

o-8. o = crs(a = 1) = / da 

Jo a 

5 m (a = l) [ a f_a8(x) 

o m (a) = / dx 

cts,o Jo x 

= /(78(a) = fa8(a) 
fa8{a « 1) a a (3.15) 



The actual data that we used in the present analysis and their corresponding references are shown in Table |III| 
We should note that some of the data in Table |III| come from surveys that overlap in the same volumes of space and 
this would in principle induce correlations of the data points. More specifically, the proper way to analyze these data 
would be to take into account their covariance matrix and calculate a total chi 2 , as was the case in the BAO data 
mentioned above. However, since this covariance matrix is not currently available 3 we perform the analysis by using 
the data "as is" , ie without a covariance matrix. 

Finally, by applying the GA on the data we will directly determine the parameter fa8(z) = GA^(z). In Fig. [5] we 
show the difference between the GA best-fit GA 3 (z) and the best-fit (f2 m = 0.27) ACDM model /ct8acdm(z). The 
gray shaded region represents the la error of the GA best-fit. In Fig. [6] we show the evolution of /cr8 in terms of the 
scale factor a and its la error region (gray shaded area). We also show fa8 for various values of w £ [—1-3, —0.7] 
(also indicated by the arrows) and f2 m = 0.27. The two vertical dotted lines indicate the regions where we trust our 
reconstruction the most 0.15 < z < 0.8. Notice that there is a sweet-spot at z ~ 1.9 (indicated by an arrow) due to 
degeneracies in fa8 for the w =const. model that, even if we had data at these high redshifts, limits its ability to 
exclude parts of the parameter space. Thus, this justifies and enhances the importance of our bias-free reconstruction 
by using GAs. Finally, clearly the best region to look for deviations from w — const, models and perhaps detecting 
modifications of gravity is clearly at z > 0.8, since as it can bee seen in Fig. [6] at this redshift region the difference 
between the different w = const, models and the extrapolation from the GA best-fit becomes maximal. 

In Fig. m we show the normalized growth factor 8(a)/5(l). The dashed line corresponds to a flat matter dominated 
model (5(a) = a, the solid black line to a ACDM model with Q m = 0.27, the dotted line to an open model OCDM 
with fi m = 0.27 and the black dot-dashed to the best-fit GA model. The gray shaded region corresponds to the la 



error region and the reason that it looks rather squashed is that, as seen in Eq. (3.151, we had to normalize 8(a) to 
unity at a = 1. 



3 The full covariance matrix between growth rate data-points at different redshifts will in principle be available in future surveys like DES. 



13 



z 


fcr& bs 


Refs. 


0.17 


0.510 ±0.060 


1201 

1 1 


0.35 


0.440 ± 0.050 


38 


0.77 


0.490 ± 0.180 


1371 


0.25 


0.351 ± 0.058 


ED 

I— ^ — j 


0.37 


0.460 ± 0.038 


i4n 


0.22 


0.420 ± 0.070 


1431 


0.41 


0.450 ± 0.040 


[431 

1 1 


0.60 


0.430 ± 0.040 


1 1 


0.78 


0.380 ± 0.040 


143] 


0.30 


0.407 ±0.055 


44 


0.40 


0.419 ± 0.041 


44. 


0.50 


0.427 ± 0.043 


BU 


0.60 


0.433 ±0.067 


m 


0.57 


0.451 ± 0.025 





TABLE III: The growth rate data and their corresponding references. 
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z 

FIG. 5: The difference between the GA best-fit f(z) and the ACDM model for fi m = 0.27. The gray shaded region represents 
the lcr error of the GA best-fit. 

IV. METHOD AND RESULTS 

In the next subsections we will apply our GA method to various models, from generic FRW models of Dark Energy 
to Large Voids of Gigaparsec size, with approximate LTB metric, and modified gravity theories like f(R). We will 
write the various observables as specific functions of redshift, and will then constrain their functional form with the 
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FIG. 6: The evolution of /<t8 in terms of the scale factor a and its la error region (gray shaded area). We also show fa8 for 
various values of w G [—1.3,-0.7] (also indicated by the arrows) and f2 m = 0.27. The two vertical dotted lines indicate the 
regions where we trust our reconstruction the most 0.15 < z < 0.8. Notice that there is a sweet-spot at z ~ 1.9 (indicated by 
an arrow) due to degeneracies in fa8 for the w =const. model that, even if we had data at these high redshifts, limits its ability 
to exclude parts of the parameter space. Thus, this justifies and enhances the importance of our bias-free reconstruction by 
using GAs. Finally, clearly the best region to look for deviations from w = const, models and perhaps detecting modifications 
of gravity is clearly at z > 0.8, since as it can bee seen in the plot, at this redshift region the difference between the different 
w = const, models and the extrapolation from the GA best-fit becomes maximal. 



GA algorithm that best fits the observational data. 



A. Generic FRW/DE models 



Taking into account all the ways we can extract information from the data, as discussed in the previous section, 
we now perform the exercise to reconstruct an FRW Dark Energy model. In particular we are interested in the 
reconstruction of the Dark Energy equation of state w(z) and the deceleration parameter q{z). The equation of state 
in given by 

, x 1 rfln(ff 2 (z)-g 2 Sl m (l + z) 3 ) (AU 
W ^ = - 1 + dln(l + z)3 (41) 



which can also be written as 



dln(l + z) 5 
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a 

FIG. 7: The normalized to unity growth factor 5(a). The dashed line corresponds to a flat matter dominated model 8(a) = a, 
the solid black line to a ACDM model with Q m = 0.27, the dotted line to an open model OCDM with f2 m = 0.27 and the black 
dot-dashed to the best-fit GA model. The gray shaded region corresponds to the la error region and the reason that it looks 
rather squashed is that, as explained in the text, we had to normalize S(a) to unity at a — 1. 



while the deceleration parameter is given by 

, . a d\nH(z) 

ff <*> s -^ = - 1+ SMiri (43) 

In order to accommodate other models like the LTB we will also follow a more geometrical approach for q(z) and 
define the deceleration parameter as 



q(z) = -1 + u"D„ ( |) = ^ ( - w^w^ + R^vTu" ) , (4.4) 



where 6 = 8 M p = D^u^ is the trace of the congruence of geodesies, where cr^w^J) is the shear (vorticity) tensor 
of the same congruence, and we have used the Raychaudhuri equation in the second equality. In a FRW space-time, 
both shear and vorticity vanish, but in a LTB model the background shear is non-zero and the first term in the RHS 
of (4.4) becomes 2(H T - H L ) 2 /(2H T + H L ) 2 . 

We will also consider the Om statistic which is defined as [50 

Om(z) = (ijy^ (4.5) 

which is constant and equal to f2 m only when H(z) corresponds to a ACDM model 

H(zf/Hl = On(l + zf + \- fi m (4.6) 

but evolves with z in a characteristic manner for different DE models, see Ref. |50| for more details. Therefore, it can 
used as a diagnostic if H(z) is determined independently, as it was shown in Ref. |20j . 

Our reconstruction method requires an a priori known value for f2 m for w(z) but fortunately this is not the case 
for q(z), which can be directly estimated by direct differentiations from the best-fit GA (Il(z). We show plots of all 
these important functions (w(z),q(z),d L (z),Om(z)) in Figs. |8{ [9| [l0| and 11 As is easily seen in Fig. [8] the results 



are consistent with a ACDM model with f2 m = 0.27, however the reconstruction is much more precise in the case of 
q(z) than w(z). Also, the Om statistic while it is consistent with a ACDM model with O m = 0.27, suffers from large 
errors especially in small redshifts, which is in agreement with the results found in |20j . 

Therefore, out of all the different diagnostics, the deceleration parameter q(z) seems to provide the best constraints 
of the cosmology assuming our model independent approach. 
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FIG. 8: Left: The evolution of w(z) with redshift for the GA best-fit and its lcr error region (gray region), both assuming 
Qm = 0.27. Right: The evolution of the deceleration parameter q(z) for the GA best-fit with its lcr error region (gray region) 
and the best-fit (fi m = 0.27 ± 0.03) ACDM model with its lcr error region (light-gray region). 
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FIG. 9: The evolution of w(z) with redshift for different values of Q m . 

B. Inhomogeneous Universes/ LTB models 

Taking into account all the ways we can extract information from the data, as discussed in the previous section, we 
now perform the exercise to reconstruct an LTB model. For more details on the LTB models see Refs. [T4lfT7] . 

The Lemaitre-Tolman-Bondi metric describes spaces with maximally symmetric (spherical) 2-dimensional surfaces, 
and is given by 



ds 2 



-dt l 



(4.7) 



for a matter source with negligible pressure and no anisotropic stress (T^, = — t)Sq S„). The function A(t, r) 
acts as an r-dependent scale factor. It is easy to see that in this framework the rates of expansion in the longitudinal 
(r) and transverse (6*, (j>) directions are, in general, different (Ht = A/ A, Hl — A' /A'). In this context, the Einstein 
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FIG. 10: Up: The evolution of q(z) with redshift for the GA best-fit and the ACDM model for different values of Q m . Down: 
The difference of q(z) between the GA best-fit and the ACDM model for different values of Q m . 



equations can be written as an effective Friedmann equation for a fixed r: 



H 2 T (t,r) = H*(r) 



(4- 



where Aq(t) = A(to,r) can be gauged to Ao(r) = r, Ho(r) = HT(to,r) and VIm{t) is the ratio between the average 
matter density inside a sphere of radius r and the critical density at that radius, and acts as an effective r-dependcnt 
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FIG. 11: The Om statistic as a function of redshift. The solid black line corresponds to the best-fit while the gray-shaded area 
to the la error region. 



matter parameter. We can then integrate the light-cone equations, 

dt = A'(r,t) dr = y/l - k(r) 

dN A'{r,t) ' dN A'(r,t) ' 

where N = ln(l + z) is the effective number of e-folds before the present time, to obtain the redshift dependence of 
all observables. 

An important parameter in the LTB models is £(z), given by Eq. (3.6). Under the assumption that the time t e is 
calculated at a high enough redshift, eg z > 100, the LTB metric will be sufficiently close to an FRW metric and we 
can simplify £ quite significantly. More specifically, we have 

A'(roo,t e ) a te _ 1 



A{Too,t e ) _ T-qoO^ _ 1 

a 1 + z e 

r 

A(r,t e ) ps ra e 



l + z e 

A'(r,t e ) « a e = —J— (4.10) 



and £ simplifies to 



while d z becomes 



(A') 1 / 3 A 2 ? 3 



1/3 



d z {z) = (l + z) &z) #4°, = ( ±±1) Ibao Hl 



1 + Z\ 1/3 l BAO /.A 1/3 



-(«) 2 /3 



A' (4.12) 
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From Eq. (4.12) we can get an expression for A' in terms of d 

A'(z)=r{zf 



d z {z) 



BAO 



l + z 



(4.13) 



Finally, the last missing ingredient, the information about the mass radial function U,m{t), will come from the 
growth rate data. The analytical solution of the growth rate equation for the LTB models, under the assumption of 
negligible shear (below a few percent), was found in Ref. [3T] to be 



JK ' 7 K m " 2F1 [1,2, 7/2,1 -n m \z)] 



where 



i-n- 1 (z) = (i-n- M 1 {r{z))) 



r(z) 



(4.14) 



(4.15) 



and f2jvf(r) is the mass radial function of the LTB model and is the desired parameter. For ease of reference we also 
provide all the relevant equations below: 



(l + z) 2 A(z) = GA 1 (z) 



A'(z) = r{zf 



d z (z) 



1 



Jbao 

n m 1 (z) = (i-n-i(r(z))) 



l + z 

A{z) 
r{z) 



Ho(r) 



(1 + z)r'{z)A' = (1 + r 2 H Q (r) 2 (l - il M (r))) 
1 fiitf(r) 



1/2 



l-n M (r) (l-fi M (r)) 3 / ; 



sinh 



(4.16) 
(4.17) 

(4.18) 

(4.19) 
(4.20) 



where Eq. (4.16) follows from the definition of the luminosity distance in an LTB model, Eq. (4.17) follows from Eq. 
(|4.13|) and Eq. (|4.18[) follows from Eq. (|4.15|). Finally, Eq. (|4. 19h follows from the redshift equation Eq. (|4.9|), see 



also |14j. while Eq. (4.20) is valid only for constrained LTB models where the Big Bang is homogeneous. All of the 
above equations can be combined to give a non-linear differential equation of the form 4 



r '( z ) = f( r ,z) 



(4.21) 



which can be solved numerically and provide the solution r — r(z). Finally, we can now calculate the mass radial 
function Cl^^r) with the following procedure: 



1. First, we calculate fcr8(z) = GA 3 (z) by applying the GA to the growth rate data of Table III 

2. Then we use Eq. ( |3.15 1 to get the growth rate f(z) and and the growth factor S m (z). 

3. With knowledge of f(z) we can now calculate the matter density parameter fl m (z) by numerically inverting 

Eq. u.uy 



4. And finally we can estimate the mass radial function fijw(r) by using Eq. (4.18), where r(z) is known from 
Eq. (|4.21|) and A(z) from Eq. (|4 16[). 



The results of this exercise are shown in Figs. 12 , 13 and |14| Clearly our reconstruction of ^Im{t) is in disagreement 
with the minimum-x 2 LTB models found in Ref. |29j . which may be interpreted as an indication that the models 
for the mass radial function Qm(t) mostly used in the literature today are clearly disfavored when a more bias-free 
method is used. 



4 We do not explicitly write down the RHS of this equation as it is quite complicated, but the interested reader may easily reproduce it 
in a program like Mathcmatica by using the relevant equations. 
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FIG. 13: as a function of r for three different cases: The best-fit from the Genetic Algorithm reconstruction (solid 

black line), a constrained CGBH (red line) and an open constrained CGBH model (green line). The solid curves correspond 
to the values of the minimum-x 2 model, while the dotted and dashed curves correspond to the marginalized and BAO+CMB 
combinations for the open and flat cases respectively. The parameter values for all these plots can be found in Table 3 of Ref. 
|29| . The two vertical dashed lines correspond to z = 0.1 and z = 0.78 and indicate the range where our reconstruction is more 
accurate (see the text for a discussion). 
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FIG. 14: -Ho(r) as a function of r for three different cases: The best-fit from the Genetic Algorithm reconstruction (solid black 
line), a constrained CGBH (red line) and an open constrained CGBH model (green line). The solid curves correspond to the 
values at the minimum while the dotted and dashed curves correspond to the marginalized and BAO+CMB combinations 
for the two cases respectively. The values for all these plots can be found in Table 3 of Ref. |29| . The two vertical dashed 
lines correspond to z = 0.1 and z — 0.78 and indicate the range where our reconstruction is more accurate (see the text for a 
discussion). 



C. Modified gravity theories 



In many modified gravity theories the growth factor of matter perturbations evolves according to Eq. (3.11) but 
with a Gat that is no longer constant and instead is time and scale dependent [ 



. The reason for this is that there 
is now a propagating scalar degree of freedom which at the linear level can be considered as a modification of Newton's 
law of universal attraction ~ 1/r 2 or equivalently a modification of Newton's constant Gjv- More specifically, it can be 
shown that on sub-horizon scales, ie when k 2 3> a 2 H 2 where k is the wave-number of the modes of the perturbations 
in Fourier space, Eq. (3.11) can be written as 



5" (a) + 



g'(o) 
H(a) 



S'(a) 



3 fi m Geff(a) 
2a 5 H{a) 2 /H 2 



2 8(a) = 



(4.22) 



where primes denote differentiation with respect to the scale factor, H(a) = - is the Hubble parameter and we assume 
the initial conditions 8(0) = and 8'(0) = 1 for the growing mode. When G e s(a) = 1 we get GR as a subcase, while 
in general in modified gravity theories G e g can be time and scale dependent, ie G e s = G e g(a,k). For example, in 
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f(R) theories we have that [37] 



. , 1 l+4^m 

G oS (a,k) = — ° (4.23) 
birr 1 -f 3\m 

m = (4.24) 
F - /,* = | (4.25) 



which reduces to GR only when /(i?) = R — 2A. 



Now, by using Eq. (4.22) and since we know both H(z), from the best-fit of the Snla data GAi, and 5(a) from the 
growth rate data as mentioned in the previous section, we can get a constraint on f2 m G c ff or on G e e by assuming 
some value for O m and we show these important parameters in Fig. |15| As it can be seen, there is a 3<r deviation 
from unity at z > 0.3, which if it persists after using future high quality data could hint towards either the existence 
of deviations from General Relativity, eg f(R) theories, or the presence of background shear. 



Also, we will use the diagnostic 



where Q m is the value of the matter density as measured from other independent observations and 



i ^ 
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m,GR Jo 



= / dx g(x) (4.27) 



/cr8(l) Jo 2//ct8(1) 



Notice that Eq. (4.271 does not include the Hubble parameter H(z) or any assumption for some dark energy model at 
all. Therefore, the calculation for I can be carried out by using only the growth rate data. Therefore, if we assume 
that the value of f2 m is independently and accurately determined by other observations, then any deviation of the 
quantity from zero, clearly and uniquely identifies an evolving G e g and consequently modified gravity theories (for 
more details see Ref. [51]). Finally, notice that the diagnostic depends solely on directly measurable quantities 
and not on any model. So, the value for the parameter that we found is = — 0.22^q'^, which hints towards an 
evolving G e g in accordance with the findings of Ref. [51] . The improvement to the constraint compared to the one 
found in |51] by one of the authors, is due to two reasons: first we are now using many more growth rate data and 
second, we now no longer demand a smooth and slowly evolving G c g , as the GA can accommodate any kind of model. 
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FIG. 16: The parameter rj( z) fr om the Snla and BAO data (dark gray area). We also show the best-fit from Ref. for 

at the 68.3% confidence (light gray area). The dotted line stands for 



the parametrization of Eq. (|5.3|) with e 

n(z) = i. 
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V. TESTS OF THE ETHERINGTON RELATION 

In this section we present results of the analysis for the Etherington relation di(z) = (1 + z) 2 <Ia{z) and the 
reconstruction of the parameter 



n 



d L (z) 



{l + zYd A {z) 



(5.1) 



The Etherington relation holds for all metric theories of gravity where photons travel along unique null geodesies and 
it depends only on the photon number conservation. Therefore, in metric theories we have that rj = 1 at all rcdshifts. 
By using the definition of the d z parameter we can rewrite rj(z) as 



_ Dl{z) ( d z (z) 



3/2 



Z 1 ' 2 

l + z 



(5.2) 



where Dl(z) and H(z) are calculated from the best-fit of the Snla data GA\ and d z (z) from the BAO data. In Fig. 
16 we show the result for the reconstruction of the parameter r/(z) for S! m = 0.27 and we also compare our results 
with the ones found in Ref. [36] , where the authors used a rather restrictive parametrization 



r,{z) = {l + zY 



(5.3) 



The best-fit they found was e = — O.Olj^QQg at the 95% confidence. The result of our reconstruction is clearly 
compatible with unity (77 = 1) despite a 3a bump at z ~ 0.5. 



VI. CONCLUSIONS 



Understanding the nature of Dark Energy is one of the most fundamental issues in Cosmology today. Several 
probes have been designed to search for deviations from a simple cosmological constant, A, and a host of observations 
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will be performed in the near future with various surveys like DES (55), PAU [S3] and, in the future, LSST [53] and 
EUCLID [55 , in order to distinguish DE from A. It is thus desirable to have the appropriate statistical tools to 
address that question in a model independent and unbiased way. 

We used Genetic Algorithms as a novel approach to extract information from several cosmological probes, such 
as type la supernovae (Snla), the Baryon Acoustic Oscillations (BAO) and the growth rate of matter perturbations. 
This information consists of a model independent and bias-free reconstruction of the various scales and distances that 
characterize the data, like the luminosity di,(z) and the angular diameter distance (Ia(z) in the Snla and BAO data 
respectively or the dependence with redshift of the matter density f2 m (a) in the growth rate data. This cosmological 
information can then be used to reconstruct the expansion history of the Universe and the resulting DE equation 
of state w(z) in the context of FRW models, or the mass radial function VLm(t) in LTB models. In this way, 
the reconstruction was completely independent of our prior bias towards specific DE equations of state w(z) or mass 
radial functions f2j\/(r). Furthermore, we used this method to test the Etherington relation, ie the well-known relation 
between the luminosity and the angular diameter distance, i] = ( 1+ ^2^( z ) i which is equal to 1 in metric theories of 
gravity. 

More specifically, in the context of FRW/DE models, we compared many different diagnostics like the equation 
of state w(z), the deceleration parameter q(z) and the Om statistic and we found that the best in constraining the 
evolution and properties of DE is the deceleration parameter q(z) as it requires no prior assumptions or a priori set 
parameters like w(z) does and has the smallest errors especially at small redshifts. Our result was in agreement with 
the ACDM model for il m = 0.27. In the case of the LTB models, we reconstructed the mass radial function flM(r) 
and the Hubble parameter H (r) and we found that they are in some tension with the commonly used profiles found 
in the literature. This of course does not mean that the LTB model is ruled out, but instead it means that these 
specific profiles are not a good fit to the latest data. 

Regarding the modified gravity theories, we reconstructed the important parameter G fj, which is a smoking gun 
signature for deviations from GR if it is different from unity at any redshift, and we found a 3er bump at intermediate 
redshifts, ie at z ~ 0.5. We found a similar deviation in the same redshift range for the Etherington relation, however 
since the statistical significance for both of these was low we cannot draw any strong conclusions. One explanation, 
at least for the anomaly in G e ff, could be the presence of a small background shear 5 , however taking this properly 
into account is beyond the scope of the current work, so we leave it for the future. 

Finally, we presented a novel way to estimate analytically the errors on the reconstructed by the Genetic Algorithm 
quantities by calculating a path integral over all possible functions that may contribute to the likelihood. We did 
this for both cases where the data are correlated and uncorrelated with each other, the latter obviously being just a 
special case of the former. In order to assess how rigorous the path integral method is, we also considered an explicit 
example and compared it numerically with a Bootstrap Monte Carlo and the Fisher Matrix approach. We found that 
the agreement between all three methods is remarkably good, something which lends support to our newly proposed 
method for the error estimation in the GA paradigm. Furthermore, we should note that since only the Bootstrap 
Monte Carlo is applicable to the case of the GAs (as the best-fit has no parameters over which one may differentiate, 
like one would do with the Fisher Matrix), our method has the obvious advantage that it is much faster that the 
Bootstrap, less computationally expensive and equally reliable. 
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